*	See Gibbs estimation of microstructure models: Teaching notes
	November, 2006
	Joel Hasbrouck
____________________________________________________________________________________________________

	RollGibbsBeta 01.sas

____________________________________________________________________________________________________;

*	Modify the following statement to set current directory to location of program;
x 'cd c:\Active\Mkt Microstructure Book\Programs Fall 2006';
options nodate nocenter nonumber ps=70 ls=120  sasautos=('./Sas Macros','.') mprint; 
libname SASData '.\Datasets';
libname this '.';

/*ods pdf file='./test.pdf' style=serifPrinter;*/

%let Infinity=1e30;
%let eps=1e-30;

*___________________________________________________________________________________________________

	Generate some price data
____________________________________________________________________________________________________;
%let nObs=100;
%let sdu = .01;
%let c = 0.01;
%let beta = 1.1;
%let sdm = 0.01;
data sim;
	m = 0;
	pm = 0;	*	pm is the market index;
	do t=1 to &nObs;
		m = m + &sdu*normal(345699);
		pm = pm + &sdm*normal(12312);
		q = 1-2*(uniform(246)>.5);
		p = m + &beta*pm + &c*q;
		keep q p pm;
		output;
	end;
	run;

%let nSweeps=10000;


*___________________________________________________________________________________________________

	Simulate parameters and q's
____________________________________________________________________________________________________;
%RollGibbsBeta(sim, nSweeps=&nSweeps, cStart=&c, varuStart=%sysevalf(&sdu*&sdu),
	betaStart=&beta, qDraw=1);

data olsIn;
	set sim;
	r = dif(p);
	rm = dif(pm);
	run;
proc means data=olsin;
proc reg data=olsIn;
	model r = rm / noint;
	quit;
	run;

*___________________________________________________________________________________________________

	Analyze the simulated estimates
____________________________________________________________________________________________________;
%let pDrop=0.2;							*	Proportion of sweeps to drop (burn in);
%let nDrop = %sysevalf(&pDrop*&nSweeps);	*	Number of sweeps to drop;
data parmOut2;
	set parmOut;
	if sweep>&nDrop then output;
	run;
proc transpose data=parmOut2 out=parmOut2T (rename=(col1=value));
	by sweep;
	run;
proc means data=parmOut2 fw=6;
	run;
proc sort data=parmOut2T;
	by _name_ sweep;
	run;
proc kde data=parmOut2T;
	univar value / levels out=myOut;
	by _name_;
	run; 
goptions ftext=complex htext=2 colors=(black) hby=0;
symbol1 v=none i=sm w=2;
proc gplot data=myOut;
	plot density*value=1;
	title "Posterior for #byval1";
	by _name_;
	run;
	quit;

/*ods pdf close;*/
